Load Packages

library(here)
library(mclust)
library(raster)
library(rgdal)
library(sp)
library(RColorBrewer)
library(PCAtools)
library(sf)
library(tmap)
library(terra)
library(rgeos)
library(spData)
library(knitr)

Import Data

#import raw data
dat_raw <- read.csv(here("data/CommunityData-raw-2015.csv"))
#select columns of interest, removed co2_hhs
dat  <- dat_raw[,c(4:12,14:21)]
#convert to dataframe
dat <- as.data.frame(unclass(dat))
#set rownames to city names
rownames(dat)=dat_raw$ME
# remove any records that have NAs
dat = na.omit(dat)
# convert to z-scores
dat <- as.matrix(scale(dat, center = TRUE, scale = TRUE))

Run Cluster GMM

Assumptions: - MSAs form clusters characterized by a multivariate distribution - Model forms: shape, volume, orientation - GMM fits a series of models with different forms and numbers of clusters - Models with highest probability and fewest parameters selected as most optimal - Based on Bayesian Information Criterion (BIC)

#run mclust on 1 to 20 clusters
dat_mc <- Mclust(dat, G=c(1:20))
#print summary of model fit
summary(dat_mc)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VVE (ellipsoidal, equal orientation) model with 10 components: 
## 
##  log-likelihood   n  df       BIC       ICL
##       -16206.34 886 485 -35704.23 -35902.49
## 
## Clustering table:
##   1   2   3   4   5   6   7   8   9  10 
##  50  59  53  62 150 129  51  66 142 124

Model comparison

#plot BIC scores of different models
plot(dat_mc, what='BIC')

#zoom in on best-fitting models
plot(dat_mc, what='BIC', ylim=c(-37000,-35000))

#model type VVE chosen as best, with similar BIC scores for 9-12 clusters
#to do additional model comparisons, generate and compare VVE models with 9-12 clusters
dat_mc9 <- Mclust(dat, G=9, modelNames="VVE")
dat_mc10 <- Mclust(dat, G=10, modelNames="VVE")
dat_mc11 <- Mclust(dat, G=11, modelNames="VVE")
dat_mc12 <- Mclust(dat, G=12, modelNames="VVE")

#extract BIC scores
BICs<-c(dat_mc9$BIC[1],dat_mc10$BIC[1],dat_mc11$BIC[1],dat_mc12$BIC[1])
#plot BIC scores
plot(c(9:12),BICs, pch=16, ylab="BIC Score",xlab="Number of Clusters", type='o')

#calculate change in BIC score, since in this method the goal to to maximize BIC
#we calculate change from highest scoring model
delta_BIC <- max(BICs) - BICs
#calculate BIC weights
w_BIC <- round(exp(-0.5*delta_BIC)/sum(exp(-0.5*delta_BIC)), digits=3)
#make table of results
results_table <- cbind.data.frame(clusters=c(9:12),BIC=BICs,delta_BIC,weight=w_BIC)
#order by delta
results_table <- results_table[order(delta_BIC),] 
rownames(results_table)<-NULL
#print table
kable(results_table[1:4], caption = "Model Comparison")
Model Comparison
clusters BIC delta_BIC weight
10 -35704.23 0.00000 1
12 -35731.46 27.23133 0
11 -35811.06 106.82569 0
9 -35818.70 114.46901 0

Extract and map cluster classifications

#extract classifications, e.g., which MSA belongs to which cluster, write to CSV
write.csv(dat_mc10$classification,file=here("results/cluster_assignment.csv"))
data<-read.csv("results/cluster_assignment.csv")

#import shapefile
msa_Boundary <-readOGR(here("data/MSA"),"tl_2015_us_cbsa") 
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\dinap\Desktop\Research\sustainable_communities\data\MSA", layer: "tl_2015_us_cbsa"
## with 929 features
## It has 12 fields
## Integer64 fields read as strings:  ALAND AWATER
#merge them
merged <- merge(msa_Boundary,data,by.x="NAME",by.y="X")
#remove any cases MSAs with no cluster assignments
merged_clean <- merged[!is.na(merged$x),]
#convert x to factor
merged_clean$Cluster <- as.factor(merged_clean$x)
#convert to sf object
cluster_sf <- st_as_sf(merged_clean)
#make us map
us_map <- tm_shape(us_states) + tm_borders()
#combine them
cluster_map <- tm_shape(cluster_sf) + tm_fill(col="Cluster", palette="Spectral")
#plot them
us_map + cluster_map

Examine how clusters score in different metrics

#new data frame
dat_clust <- as.data.frame(dat)
#add names
dat_clust$ME <- rownames(dat)
rownames(dat_clust) <- NULL
#merge in cluster assignments
dat_clust <- merge(dat_clust,data,by.x="ME",by.y="X")
#calculate cluster means for different variables
c1m <- data.frame(one=colMeans(subset(dat_clust, x=="1")[2:18]))
c2m <- data.frame(two=colMeans(subset(dat_clust, x=="2")[2:18]))
c3m <- data.frame(three=colMeans(subset(dat_clust, x=="3")[2:18]))
c4m <- data.frame(four=colMeans(subset(dat_clust, x=="4")[2:18]))
c5m <- data.frame(five=colMeans(subset(dat_clust, x=="5")[2:18]))
c6m <- data.frame(six=colMeans(subset(dat_clust, x=="6")[2:18]))
c7m <- data.frame(seven=colMeans(subset(dat_clust, x=="7")[2:18]))
c8m <- data.frame(eight=colMeans(subset(dat_clust, x=="8")[2:18]))
c9m <- data.frame(nine=colMeans(subset(dat_clust, x=="9")[2:18]))
c10m <- data.frame(ten=colMeans(subset(dat_clust, x=="10")[2:18]))
cluster_means <- cbind.data.frame(c1m,c2m,c3m,c4m,c5m,c6m,c7m,c8m, c9m, c10m)
cluster_means$variable <- row.names(cluster_means)
rownames(cluster_means) <- NULL
#boxplots of scores for clusters
par(mfrow=c(2,5), mar=c(2.5,5,2.5,3))
boxplot(subset(dat_clust, x=="1")[2:18], horizontal=T, las=1, ylim=c(-4,4), main="Cluster 1")
boxplot(subset(dat_clust, x=="2")[2:18], horizontal=T, las=1, ylim=c(-4,4), main="Cluster 2")
boxplot(subset(dat_clust, x=="3")[2:18], horizontal=T, las=1, ylim=c(-4,4), main="Cluster 3")
boxplot(subset(dat_clust, x=="4")[2:18], horizontal=T, las=1, ylim=c(-4,4), main="Cluster 4")
boxplot(subset(dat_clust, x=="5")[2:18], horizontal=T, las=1, ylim=c(-4,4), main="Cluster 5")
boxplot(subset(dat_clust, x=="6")[2:18], horizontal=T, las=1, ylim=c(-4,4), main="Cluster 6")
boxplot(subset(dat_clust, x=="7")[2:18], horizontal=T, las=1, ylim=c(-4,4), main="Cluster 7")
boxplot(subset(dat_clust, x=="8")[2:18], horizontal=T, las=1, ylim=c(-4,4), main="Cluster 8")
boxplot(subset(dat_clust, x=="9")[2:18], horizontal=T, las=1, ylim=c(-4,4), main="Cluster 9")
boxplot(subset(dat_clust, x=="10")[2:18], horizontal=T, las=1, ylim=c(-4,4), main="Cluster 10")

par(mfrow=c(1,1))
#dotplots of cluster mean values
par(mfrow=c(2,5))
dotchart(cluster_means$one, labels=cluster_means$variable,xlim=c(-2,2), main='Cluster 1', pch=16)
abline(v=0, lwd=2)
dotchart(cluster_means$two, labels=cluster_means$variable,xlim=c(-2,2), main='Cluster 2', pch=16)
abline(v=0, lwd=2)
dotchart(cluster_means$three, labels=cluster_means$variable,xlim=c(-2,2), main='Cluster 3', pch=16)
abline(v=0, lwd=2)
dotchart(cluster_means$four, labels=cluster_means$variable,xlim=c(-2,2), main='Cluster 4', pch=16)
abline(v=0, lwd=2)
dotchart(cluster_means$five, labels=cluster_means$variable,xlim=c(-2,2), main='Cluster 5', pch=16)
abline(v=0, lwd=2)
dotchart(cluster_means$six, labels=cluster_means$variable,xlim=c(-2,2), main='Cluster 6', pch=16)
abline(v=0, lwd=2)
dotchart(cluster_means$seven, labels=cluster_means$variable,xlim=c(-2,2), main='Cluster 7', pch=16)
abline(v=0, lwd=2)
dotchart(cluster_means$eight, labels=cluster_means$variable,xlim=c(-2,2), main='Cluster 8', pch=16)
abline(v=0, lwd=2)
dotchart(cluster_means$nine, labels=cluster_means$variable,xlim=c(-2,2), main='Cluster 9', pch=16)
abline(v=0, lwd=2)
dotchart(cluster_means$ten, labels=cluster_means$variable,xlim=c(-2,2), main='Cluster 10', pch=16)
abline(v=0, lwd=2)

par(mfrow=c(1,1))

Which MSAs best characterize their cluster?

Typical members of each cluster
Cluster Name
1 Goldsboro, NC
2 Sierra Vista-Douglas, AZ
3 Lexington-Fayette, KY
4 Tampa-St. Petersburg-Clearwater, FL
5 Texarka, TX-AR
6 Zanesville, OH
7 Louisville/Jefferson County, KY-IN
8 Watertown, SD
9 Pama City, FL
10 Williamsport, PA

Which cluster scores the best on the different metrics?

Metric Cluster
AQI_Good eight
Bachelor_Over_25 three
Per_Poverty eight
Gini eight
Per_Sev_Hous eight
Xstreamlengthimpaired two
Per_Avg_Land_Cov ten
poor_health_percent eight
Z_Water_Index six
Index_Black five
Index_Asian four
Index_Latino three
GHG_Percap four
UNEMPLOY eight
FOOD_INS eight
VIO_CRIME eight

Examine cluster patterns with PCA

Note: adjusting values for variables so that high=‘good’ and low=‘bad’ for all variables

dat_clust_adj <- dat_clust
dat_clust_adj$Per_Poverty <- dat_clust$Per_Poverty * -1
dat_clust_adj$Gini <- dat_clust$Gini * -1
dat_clust_adj$Per_Sev_Hous <- dat_clust$Per_Sev_Hous * -1
dat_clust_adj$Xstreamlengthimpaired <- dat_clust$Xstreamlengthimpaired * -1
dat_clust_adj$poor_health_percent <- dat_clust$poor_health_percent * -1
dat_clust_adj$Z_Water_Index <- dat_clust$Z_Water_Index * -1
dat_clust_adj$Index_Black <- dat_clust$Index_Black * -1
dat_clust_adj$Index_Asian <- dat_clust$Index_Asian * -1
dat_clust_adj$Index_Latino <- dat_clust$Index_Latino * -1

cluster1 <- data.frame(subset(dat_clust_adj, x=="1"))
cluster2 <- data.frame(subset(dat_clust_adj, x=="2"))
cluster3 <- data.frame(subset(dat_clust_adj, x=="3"))
cluster4 <- data.frame(subset(dat_clust_adj, x=="4"))
cluster5 <- data.frame(subset(dat_clust_adj, x=="5"))
cluster6 <- data.frame(subset(dat_clust_adj, x=="6"))
cluster7 <- data.frame(subset(dat_clust_adj, x=="7"))
cluster8 <- data.frame(subset(dat_clust_adj, x=="8"))
cluster9 <- data.frame(subset(dat_clust_adj, x=="9"))
cluster10 <- data.frame(subset(dat_clust_adj, x=="10"))

Perform PCA on individual clusters

#PCA cluster 1
c1_PCA <- pca(cluster1[2:18],transposed = T) #rotate table and only use numeric columns
#scree plot
screeplot(c1_PCA, title="Scree plot Cluster 1") 
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).

#elbow method of estimating important components
findElbowPoint(c1_PCA$variance) 
## PC5 
##   5
#biplot for PC1 and PC2
biplot(c1_PCA, showLoadings = T, showLoadingsNames = T, hline=0,vline=0, title='Cluster 1')
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_label_repel).

#loadings for components
plotloadings(c1_PCA, title="Cluster 1")
## -- variables retained:
## GHG_Percap, UNEMPLOY, FOOD_INS, Per_Poverty, Index_Asian, Z_Water_Index, Per_Sev_Hous, Per_Avg_Land_Cov, VIO_CRIME

#PCA cluster 2
c2_PCA <- pca(cluster2[2:18],transposed = T)
screeplot(c2_PCA)
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).

findElbowPoint(c2_PCA$variance)
## PC6 
##   6
biplot(c2_PCA, showLoadings = T, showLoadingsNames = T, hline=0,vline=0, title='Cluster 2')
## Warning: Removed 2 rows containing missing values (geom_segment).
## Warning: Removed 2 rows containing missing values (geom_label_repel).

plotloadings(c2_PCA, title="Cluster 2")
## -- variables retained:
## UNEMPLOY, Per_Poverty, Index_Asian, non_migration, Z_Water_Index, VIO_CRIME, AQI_Good, FOOD_INS, Bachelor_Over_25

#PCA cluster 3
c3_PCA <- pca(cluster3[2:18],transposed = T)
screeplot(c3_PCA)

findElbowPoint(c3_PCA$variance)
## PC4 
##   4
biplot(c3_PCA, showLoadings = T, showLoadingsNames = T, hline=0,vline=0, title='Cluster 3')
## Warning: Removed 2 rows containing missing values (geom_segment).

## Warning: Removed 2 rows containing missing values (geom_label_repel).

plotloadings(c3_PCA, title="Cluster 3")
## -- variables retained:
## Per_Poverty, non_migration, FOOD_INS, Bachelor_Over_25, Per_Avg_Land_Cov, Per_Sev_Hous, GHG_Percap, Index_Asian, VIO_CRIME, AQI_Good, Index_Latino

#PCA cluster 4
c4_PCA <- pca(cluster4[2:18],transposed = T)
screeplot(c4_PCA)
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).

findElbowPoint(c4_PCA$variance)
## PC4 
##   4
biplot(c4_PCA, showLoadings = T, showLoadingsNames = T, hline=0,vline=0, title='Cluster 4')
## Warning: Removed 2 rows containing missing values (geom_segment).
## Warning: Removed 2 rows containing missing values (geom_label_repel).

plotloadings(c4_PCA, title="Cluster 4")
## -- variables retained:
## Xstreamlengthimpaired, AQI_Good, Index_Latino, UNEMPLOY, Bachelor_Over_25, Per_Sev_Hous, Index_Black, Per_Avg_Land_Cov, Gini

#PCA cluster 5
c5_PCA <- pca(cluster5[2:18],transposed = T)
screeplot(c5_PCA)
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).

findElbowPoint(c5_PCA$variance)
## PC6 
##   6
biplot(c5_PCA, showLoadings = T, showLoadingsNames = T, hline=0,vline=0, title='Cluster 5')
## Warning: Removed 3 rows containing missing values (geom_segment).
## Warning: Removed 3 rows containing missing values (geom_label_repel).

plotloadings(c5_PCA, title="Cluster 5")
## -- variables retained:
## Per_Poverty, FOOD_INS, VIO_CRIME, Index_Black, non_migration, poor_health_percent, Index_Latino, UNEMPLOY, Index_Asian

#PCA cluster 6
c6_PCA <- pca(cluster6[2:18],transposed = T)
screeplot(c6_PCA)

findElbowPoint(c6_PCA$variance)
## PC10 
##   10
biplot(c6_PCA, showLoadings = T, showLoadingsNames = T, hline=0,vline=0, title='Cluster 6')
## Warning: Removed 2 rows containing missing values (geom_segment).
## Warning: Removed 2 rows containing missing values (geom_label_repel).

plotloadings(c6_PCA, title="Cluster 6")
## -- variables retained:
## Index_Latino, Xstreamlengthimpaired, UNEMPLOY, Gini, poor_health_percent, Index_Asian, non_migration, FOOD_INS, Index_Black, VIO_CRIME, Bachelor_Over_25, Per_Avg_Land_Cov

#PCA cluster 7
c7_PCA <- pca(cluster7[2:18],transposed = T)
screeplot(c7_PCA)
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).

findElbowPoint(c7_PCA$variance)
## PC6 
##   6
biplot(c7_PCA, showLoadings = T, showLoadingsNames = T, hline=0,vline=0, title='Cluster 7')
## Warning: Removed 2 rows containing missing values (geom_segment).
## Warning: Removed 2 rows containing missing values (geom_label_repel).

plotloadings(c7_PCA, title="Cluster 7")
## -- variables retained:
## Xstreamlengthimpaired, AQI_Good, Bachelor_Over_25, Per_Avg_Land_Cov, VIO_CRIME, Gini, poor_health_percent

#PCA cluster 8
c8_PCA <- pca(cluster8[2:18],transposed = T)
screeplot(c8_PCA)

findElbowPoint(c8_PCA$variance)
## PC9 
##   9
biplot(c8_PCA, showLoadings = T, showLoadingsNames = T, hline=0,vline=0, title='Cluster 8')
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_label_repel).

plotloadings(c8_PCA, title="Cluster 8")
## -- variables retained:
## Gini, Index_Black, GHG_Percap, non_migration, Bachelor_Over_25, VIO_CRIME, Index_Asian, UNEMPLOY, Per_Poverty, Index_Latino

#PCA cluster 9
c9_PCA <- pca(cluster9[2:18],transposed = T)
screeplot(c9_PCA)
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).

findElbowPoint(c9_PCA$variance)
## PC9 
##   9
biplot(c9_PCA, showLoadings = T, showLoadingsNames = T, hline=0,vline=0, title='Cluster 9')
## Warning: Removed 3 rows containing missing values (geom_segment).
## Warning: Removed 3 rows containing missing values (geom_label_repel).

plotloadings(c9_PCA, title="Cluster 9")
## -- variables retained:
## Index_Asian, non_migration, Per_Poverty, FOOD_INS, Bachelor_Over_25, Gini, Index_Black, Per_Avg_Land_Cov, UNEMPLOY, Per_Sev_Hous

#PCA cluster 10
c10_PCA <- pca(cluster10[2:18],transposed = T)
screeplot(c10_PCA)

findElbowPoint(c10_PCA$variance)
## PC9 
##   9
biplot(c10_PCA, showLoadings = T, showLoadingsNames = T, hline=0,vline=0, title='Cluster 10')
## Warning: Removed 2 rows containing missing values (geom_segment).
## Warning: Removed 2 rows containing missing values (geom_label_repel).

plotloadings(c10_PCA, title="Cluster 10")
## -- variables retained:
## Index_Latino, non_migration, Gini, Bachelor_Over_25, VIO_CRIME, Per_Poverty, AQI_Good, Per_Sev_Hous, Xstreamlengthimpaired, Index_Black, Per_Avg_Land_Cov

Dimension reduction

#run discriminant analysis on clusters
DR <- MclustDR(dat_mc10, lambda = 1) #setting lambda to 1 gives most separating directions
summary(DR)
## ----------------------------------------------------------------- 
## Dimension reduction for model-based clustering and classification 
## ----------------------------------------------------------------- 
## 
## Mixture model type: Mclust (VVE, 10) 
##         
## Clusters   n
##       1   50
##       2   59
##       3   53
##       4   62
##       5  150
##       6  129
##       7   51
##       8   66
##       9  142
##       10 124
## 
## Estimated basis vectors: 
##                            Dir1       Dir2       Dir3       Dir4      Dir5
## AQI_Good              -0.176001  0.0167446  0.2544551 -0.0041113 -0.199672
## Bachelor_Over_25       0.072840 -0.2897267  0.3135171 -0.5609009  0.421081
## Per_Poverty            0.022156  0.2216694  0.6106306 -0.1721047  0.485515
## Gini                  -0.018424 -0.1870715 -0.0774274 -0.0713467 -0.052922
## non_migration          0.122296  0.0306189 -0.1004147  0.4486234  0.108786
## Per_Sev_Hous          -0.096523 -0.2903113 -0.5171334 -0.2399188  0.140611
## Xstreamlengthimpaired  0.474859  0.2044703 -0.1971267  0.1582116  0.240913
## Per_Avg_Land_Cov      -0.253104 -0.3671406 -0.0276019  0.2069777 -0.338410
## poor_health_percent   -0.160339 -0.3429140 -0.0287426  0.2401105  0.097547
## Z_Water_Index         -0.751602  0.6152728 -0.1207784  0.0101883  0.270779
## Index_Black            0.012783  0.0364651 -0.1388193 -0.0924190  0.063046
## Index_Asian            0.030941 -0.1078755  0.1763719  0.0492666  0.086927
## Index_Latino           0.053899  0.0566457 -0.0075445  0.2151389 -0.028245
## GHG_Percap            -0.134751 -0.0033918  0.1659712  0.0513501  0.460770
## UNEMPLOY               0.016303 -0.0286237 -0.2089966  0.0738562  0.156282
## FOOD_INS               0.113812 -0.2359471  0.0577083  0.3324095 -0.089733
## VIO_CRIME             -0.151807 -0.0315984 -0.0236525  0.3006170  0.013413
##                            Dir6      Dir7       Dir8       Dir9
## AQI_Good               0.247547  0.396043  0.0874908 -0.0832171
## Bachelor_Over_25      -0.170414 -0.031725 -0.3655544  0.0382844
## Per_Poverty           -0.064981  0.201735 -0.6034751 -0.4234419
## Gini                  -0.079564 -0.351825  0.2765681 -0.1192998
## non_migration          0.292492 -0.165084 -0.1740190 -0.2571548
## Per_Sev_Hous           0.351465  0.364975  0.3325524 -0.2378273
## Xstreamlengthimpaired -0.438677  0.568038 -0.0946994 -0.0208607
## Per_Avg_Land_Cov      -0.188715  0.293543 -0.4090132  0.0051191
## poor_health_percent   -0.296444  0.183355  0.1365397  0.1702002
## Z_Water_Index         -0.218556  0.071558 -0.0198185  0.0456237
## Index_Black            0.105431 -0.086784 -0.1339241  0.4102571
## Index_Asian           -0.134440 -0.048924  0.0340747 -0.1862869
## Index_Latino           0.026287 -0.080645 -0.1474826 -0.1506777
## GHG_Percap             0.322789  0.095990 -0.0615401  0.1927214
## UNEMPLOY               0.388493  0.056520 -0.0461231  0.2111648
## FOOD_INS              -0.116806 -0.038561  0.1871896  0.5775028
## VIO_CRIME             -0.166452 -0.198262  0.0084581 -0.0539742
## 
##                 Dir1     Dir2     Dir3     Dir4     Dir5     Dir6     Dir7
## Eigenvalues  0.94731  0.67751  0.50022  0.43587  0.32617  0.18446  0.13214
## Cum. %      28.77575 49.35584 64.55074 77.79092 87.69873 93.30189 97.31587
##                 Dir8       Dir9
## Eigenvalues  0.06914   0.019223
## Cum. %      99.41608 100.000000
#plot eigenvalues
plot(c(1:17),DR$evalues,ylab="eigenvalues", xlab="component", type="o", pch=16)

#plot all pairs of combinations
plot(DR, what='pairs',symbols=c("1","2","3","4","5","6","7","8","9","10"),
     colors=c("red","blue","green","goldenrod2","violet","brown","black","grey30","slateblue2","seagreen3"))

#plot directions for variables
Directions <- data.frame(d1=DR$basis[,1],d2=DR$basis[,2],
                         d3=DR$basis[,3],d4=DR$basis[,4],
                         d5=DR$basis[,5],d6=DR$basis[,6],
                         d7=DR$basis[,7],d8=DR$basis[,8],
                         d9=DR$basis[,9],var=rownames(DR$basis))

par(mfrow=c(3,3))
dotchart(Directions$d1, labels=Directions$var, pch=16, main="Dir1", xlim=c(-0.6,0.6))
abline(v=0, lwd=2)
abline(v=c(0.1,-0.1),lty=2)
dotchart(Directions$d2, labels=Directions$var, pch=16, main="Dir2", xlim=c(-0.6,0.6))
abline(v=0, lwd=2)
abline(v=c(0.1,-0.1),lty=2)
dotchart(Directions$d3, labels=Directions$var, pch=16, main="Dir3", xlim=c(-0.6,0.6))
abline(v=0, lwd=2)
abline(v=c(0.1,-0.1),lty=2)
dotchart(Directions$d4, labels=Directions$var, pch=16, main="Dir4", xlim=c(-0.6,0.6))
abline(v=0, lwd=2)
abline(v=c(0.1,-0.1),lty=2)
dotchart(Directions$d5, labels=Directions$var, pch=16, main="Dir5", xlim=c(-0.6,0.6))
abline(v=0, lwd=2)
abline(v=c(0.1,-0.1),lty=2)
dotchart(Directions$d6, labels=Directions$var, pch=16, main="Dir6", xlim=c(-0.6,0.6))
abline(v=0, lwd=2)
abline(v=c(0.1,-0.1),lty=2)
dotchart(Directions$d7, labels=Directions$var, pch=16, main="Dir7", xlim=c(-0.6,0.6))
abline(v=0, lwd=2)
abline(v=c(0.1,-0.1),lty=2)
dotchart(Directions$d8, labels=Directions$var, pch=16, main="Dir8", xlim=c(-0.6,0.6))
abline(v=0, lwd=2)
abline(v=c(0.1,-0.1),lty=2)
dotchart(Directions$d9, labels=Directions$var, pch=16, main="Dir9", xlim=c(-0.6,0.6))
abline(v=0, lwd=2)
abline(v=c(0.1,-0.1),lty=2)

par(mfrow=c(1,1))